%% Description

% This script enables the user to simulate the impact of saturation on the
% steady-state and transient performance of current transformers.

%% Initialization

% Prepare the environment.

clear variables ;
close all ;
clc ;

% Define CT parameters, including its ratio, burden, internal resistance
% and saturation curve. 

f = 50 ; 
ratio = 1500 ;
Rb = 5 ;          
Rct = 3 ;       
R = Rb + Rct ;
data = [0.0010000176688569297, 10.085170417662166
    0.0019613576030532783, 26.74582417536411
    0.002966506077935653, 47.61199621748684
    0.003959045986466353, 69.73683720084738
    0.004939782146902115, 95.44258054069448
    0.005987749989810607, 122.05516704525463
    0.006984041542808529, 148.34425035006987
    0.007914283481852741, 172.81007583177555
    0.008882796724630711, 201.310952318753
    0.010260956472246778, 234.51236467793623
    0.014925127432752375, 302.45654261610525
    0.019906845942547487, 332.03100290262404
    0.029785858972400422, 334.8589248220246
    0.039722186777592726, 343.4880260518137
    0.04953066466392643, 346.4135279158924
    0.059436902629733786, 352.33949372464696
    0.06996763206110256, 355.34038390858683
    0.08002588939357885, 355.34038390858717
    0.08893475567007735, 358.3668327995606
    0.10171964822652914, 358.36683279956026
    0.1152336158789439, 361.419058082141] ;
i_o_data = data(:,1)'*sqrt(2) ;
U_data  = data(:,2)' ;
lambda_data = U_data*sqrt(2)/(2*pi*f) ;

% Define the shape of primary current, including its peak value, voltage
% phase angle at the onset of the fault, as well as the network time
% constant.

Ipeak = 10e3*sqrt(2) ;  
theta = 60*pi/180 ; 
tau = 25e-3 ;      

% Define the primary current function.

i_pp_fun = @(t) Ipeak/ratio*(sin(2*pi*f*t - theta) + sin(theta)*exp(-t/tau)) ;

% Define the inverse saturation function.

i_o_fun = @(lam) interp1(lambda_data,i_o_data,abs(lam),'pchip','extrap')*sign(lam) ;

%% ODE integration

% Define the differential equation of the current transformer:
% d(lambda)/dt = R*(i_pp - i_o)

ode_fun = @(t,lam) R*(i_pp_fun(t) - i_o_fun(lam)) ;

% Define the solver settings.

tspan = [0 0.1] ;
options = odeset('RelTol',1e-6) ;

% Solve the ODE.

[t,lambda] = ode45(ode_fun,tspan,0.7*max(lambda_data),options) ;

% Calculate currents.

i_pp = i_pp_fun(t) ;
i_o = arrayfun(i_o_fun,lambda) ;
i_s = i_pp - i_o ;

%% RMS evaluation using a Fourier filter

% Define the period and create uniform time vector.

T = 1/f;
N = 50 ; 
t_uniform = 0 : 1/(N*f) : t(end) ;
dt = t_uniform(2) - t_uniform(1) ;

% Resample the ideal and real secondary currents.

i_pp_uniform = interp1(t,i_pp,t_uniform) ;
i_s_uniform = interp1(t,i_s,t_uniform) ;

% Define the quadrature reference signals.

cos_ref = cos(2*pi*f*t_uniform) ;
sin_ref = sin(2*pi*f*t_uniform) ;

% Extract the fundamental component using moving average of quadrature
% products.

i_pp_cos = filter(ones(1,N)/N, 1, 2*i_pp_uniform.*cos_ref) ;
i_pp_sin = filter(ones(1,N)/N, 1, 2*i_pp_uniform.*sin_ref) ;
i_pp_rms = sqrt(i_pp_cos.^2 + i_pp_sin.^2)/sqrt(2) ;

i_s_cos = filter(ones(1,N)/N,1,2*i_s_uniform.*cos_ref) ;
i_s_sin = filter(ones(1,N)/N,1,2*i_s_uniform.*sin_ref) ;
i_s_rms = sqrt(i_s_cos.^2 + i_s_sin.^2)/sqrt(2) ;

%% Plot results

% Plot the current waveforms.

figure ;
plot(t*1000, i_pp, t*1000, i_s, t*1000, i_o,'--','LineWidth',2) ;
legend('Idealna sekundarna struja', 'Stvarna sekundarna struja','Struja magnećenja','Location','southeast') ;
xlabel('Vrijeme [ms]','FontWeight','bold') ; 
ylabel('Struja [A]','FontWeight','bold') ;
grid on ;
set(gcf,'color','w') ; 
set(gca,'FontName','LM Roman 10','FontSize',12) ; 

% Plot the RMS current used by protection relays.

figure;
plot(t_uniform*1000, i_pp_rms,t_uniform*1000, i_s_rms,'--','LineWidth',2) ;
legend('Idealna efektivna vrijednost', 'Stvarna efektivna vrijednost') ;
xlabel('Vrijeme [ms]','FontWeight','bold') ; 
ylabel('Efektivna vrijednost struje [A]','FontWeight','bold') ;
grid on ;
set(gcf,'color','w') ; 
set(gca,'FontName','LM Roman 10','FontSize',12) ; 

% Plot the flux linkage.

figure ;
plot(t*1000, lambda, 'LineWidth', 1.5) ;
xlabel('Time [ms]') ;
ylabel('Flux linkage [Wb-turns]') ;
grid on ;
set(gcf,'color','w') ;

%% Plot the current waveforms with RMS values.

figure ;
p1 = plot(t*1000,i_pp,'LineWidth',1.5) ;
hold on ;
p2 = plot(t*1000,i_s,'LineWidth',1.5) ;
plot(t_uniform*1000, i_pp_rms,'LineWidth',2,'Color',p1.Color,'LineStyle','--') ;
plot(t_uniform*1000, i_s_rms,'LineWidth',2,'Color',p2.Color,'LineStyle','--') ;
legend('Idealna sekundarna struja', 'Stvarna sekundarna struja',...
'Idealna efektivna vrijednost', 'Stvarna efektivna vrijednost',...
'Location','southwest') ;
xlabel('Vrijeme [ms]','FontWeight','bold') ;
ylabel('Struja [A]','FontWeight','bold') ;
grid on ;
set(gcf,'color','w') ;
set(gca,'FontName','LM Roman 10','FontSize',12) ;

